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The topological and dynamical features of small scales are studied in the context of decaying mag- 
netohydrodynamic turbulent flows using direct numerical simulations. Joint probability density 
functions (PDFs) of the invariants of gradient quantities related to the velocity and the magnetic 
fields demonstrate that structures and dynamics at the time of maximum dissipation depend on 
^y^ the large scale initial conditions. This is evident in particular from the fact that each fiow has 

^H a different shape for the joint PDF of the invariants of the velocity gradient in contrast to the 

^^ universal teardrop shape of hydrodynamic turbulence. The general picture that emerges from the 

C^ analysis of the invariants is that regions of high vorticity are correlated with regions of high strain 

;-H rate S also in contrast to hydrodynamic turbulent fiows. Magnetic strain dominated regions are 

Qh also well correlated with region of high current density j. Viscous dissipation (ex S'^) as well as 

Ohmic dissipation (ex j'^) reside in regions where strain and rotation are locally almost in balance. 
The structures related to the velocity gradient possess different characteristics than those associated 
with the magnetic field gradient with the latter being locally more quasi-two dimensional. 
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J3 I. INTRODUCTION 

q 

C/^ Numerous studies have shown that the key assumptions of small-scale universality, isotropy and locality of 

^^ interactions are in question in various contexts of magnetohydrodynamic (MHD) turbulence^ ^. Therefore, 

c/^ the validity of the classical phenomenology of Kolmogorov (K41)^, which provides to a good approximation 

^^ the power law of the energy spectrum in hydrodynamic turbulence (besides intermittency corrections), is 

^^ questionable in MHD turbulence, where several debatable phenomenological theories exist^ ^^. In summary 

\_j the power law scaling exponents obtained in the various MHD turbulence phenomenologies based on weak 

and strong turbulence arguments both for isotropic and anisotropic fields are —5/3, —3/2 and —2. 
^"^ Numerical simulations to date are unable to provide a definitive answer to the scaling of the energy 

J^ spectrum in MHD turbulence^^^. Recently, large resolution simulations by Lee et al.^^ (using a code that 

enforces the symmetries of the Taylor- Green vortex to achieve higher resolution) demonstrated different 
\^ scaling of the total energy spectra for different initial conditions and thereby suggested that freely decaying 

^^ MHD turbulent fiows are non-universal. 

I* This lack of the detailed knowledge of the energy spectrum in MHD turbulence has many implications 

^^ because to predict for example heating rates in solar and space physics, the energy dissipation rate is 

^y^ required, which is dependent on the slope of the energy spectrum. This is also the reason why subgrid scale 

i-H models, required for numerical modelling in astrophysics and geophysics, are less developed in MHD. 

h^ On the other hand, several universal small scale features have been observed in a variety of hydrodynamic 

. i-H turbulent fiows since the seminal works by Perry, Chong and CantwelP^^^ on the analysis of the velocity 

^^ gradient tensor invariants. One of the most important universal results is the well known teardrop shape^^^^ 

H of the joint probability density function (PDF) of the invariants of the velocity gradient tensor, which 

describes the topology and dynamics of small scales in hydrodynamic turbulence. Other important directions 

of research in hydrodynamic turbule nce in volving the study of such invariants has been the topological 

classification of the coherent structures^^^ and the use of the invariants in subgrid-scale modelling^^. 

Therefore, due to the limited information that can be extracted just from the slopes of the energy spectra 
and the fact that small scale universality is one of the key assumptions in inertial range phenomenologies, 
we try to gain insight in this paper on the universal or/and non-universal features of MHD turbulence by 
studying the structures and dynamics of small scales through joint PDFs of the invariants of the velocity 
gradient, magnetic field gradient and related tensors. Through this analysis we also attempt to provide a 
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classification of the structures in MHD turbulence. This investigation was carried out using Direct Numer- 
ical Simulation (DNS) data of incompressible, homogeneous, decaying MHD turbulence with no imposed 
symmetries and no magnetic flux either in or out of our periodic boxes. 

The paper is organised as follows. The numerical method, the initial conditions and the parameters of our 



DNS of decaying MHD turbulent flows are provided in sectionjilj In section [Hi] we present the energy spectra 
of our flows. Before presenting our results we give an outline for the classification of fluid flow topology 
in section [TV| Then, in sections \V\ and [Vl| we unravel our joint PDF analysis for the invariants of gradient 
quantities related to the velocity and magnetic field, respectively, delineating the structure and dynamics of 
the examined MHD flows. At the end, in section [VTH we summarise our results. 

II. DNS OF DECAYING MHD TURBULENCE 

A. Governing equations 8i numerical method 

We consider the three-dimensional, incompressible MHD equations of fluid velocity u and magnetic in- 
duction b to be 

dfU - uAu = (n X u;) + {j x b) - VP (1) 

dtb - K.Ab = V X (n X 6) (2) 

V • n = V • 6 = (3) 



with u the kinematic viscosity, tz the magnetic diffusivity, u; = V x n the vorticity, j = ^ x b the current 
density of the magnetic field and P = p/p -^ \u^ the fluid pressure, composed by the plasma pressure p, 
the constant mass density p and the hydrodynamic pressure \v?' - Note that magnetic induction has units 
of Alfven velocity, i.e. bj ^pp^^ where /io = (^cr)~^ is the permeability of free space with a the electrical 
conductivity. In ideal MHD, where i^ = z^: = 0, the total energy Et = ^(|np + |6p) = Eu-\-E}), the magnetic 
helicity H^j = {u - b) and the cross helicity He = {a - b) are conserved, where the angle brackets (.) in this 
study denote spatial averages. Here, a is the magnetic potential, which is defined as a = — A~^(V x 6), 
since one can define 6 = V x a with V • a = 0. 

Our numerical method is pseudo-spectral^^, where each component of u and b is represented as truncated 
Galerkin expansions in terms of the Fourier basis. The non-linear terms are initially computed in physical 
space and then transformed to spectral space using fast Fourier transforms^^. Aliasing errors are removed 
using the 2/3 dealiasing rule, i.e. wavenumbers k G [1, A^/3], where N is the number of grid points in each 
Cartesian coordinate of our periodic box with period 27r. The non-linear terms along with the pressure 
term are computed in such a way that u and b are projected on to a divergence- free space so that Eqs. 
PI are satisfied. The temporal integration of Eqs. M and Q is performed using a second-order Runge- 
Kutta method. The code is parallelised using message passing interface (MPI) with one-dimensional domain 
decomposition^^. 

B. Initial conditions & numerical parameters 

The initial conditions that we consider in this study are the three different cases studied iiP^. In particular, 
the initial velocity field is the Taylor-Green (TG) vortex^^ defined as 

utg {^) = '^0 (sin X cos y cos 2:, — cos x sin y cos 2:, 0) (4) 

and the initial conditions of the magnetic field are generalisation of the TG vortex symmetries. In detail, 
the insulating case (run "I" hereafter) is 

6/ {x) = 6q (cos X sin y sin z, sin x cos y sin z ^ —2 sin x sin y cos z) (5) 

where ji = ^ x bj is parallel to the faces of a subvolume [0, tt]^, which can thereby be considered as electrical 
insulators. Note that in this case the magnetic field 6/ = — (6o/iio)V x Utg and the magnetic as well as 
cross helicity are globally restricted to vanish for all times due to the TG symmetries. The conducting case 
(run "C" hereafter) takes the following form 

be {x) = 6q (sin 2x cos 2y cos 2z, cos 2x sin 2y cos 2z, —2 cos 2x cos 2y sin 2z) (6) 



with jc = ^ ^ be perpendicular to the faces of a subvolume [0,7r]^, which can consequently be considered 
as electrically conductive. In this configuration, Hjj = for all times but He ^ although negligible (i.e. 
Hci/Ef < 0.04 at its maximum over time, where ^ is a typical length scale). The final case that is considered 
by Lee et al.'^is an alternative (run "A" hereafter) to the insulating initial conditions above (see Eq. ([s])), 
namely 

bA{x) = b^ (cos 2x sin 2y sin 2z, — sin 2x cos 2y sin 2z, 0) (7) 

for which again Hij = He = for all times, at least up to the peak of dissipation. 

The above TG fields exhibit several intrinsic symmetries within a cubic box of size [0, 27r]^, where periodic 
boundary conditions are applied. These are mirror ( ant i) symmetries about the planes x = 0, x = tt, ?/ = 0, 
y = TT, z = and z = tt as well as rotational (ant i) symmetries of angle Nir about the axes (x, ?/, z) = (|, 7/, |) 
and (x, 1^, ^) and of angle N7t/2 about the axis (|^, |^, z) for A/' G Z. The above mentioned planes that possess 
mirror symmetries form the insulating and conducting walls of [0,7r]^ sub-boxes, also called impermeable 
boxes^, for the corresponding initial conditions. 

It is important to mention that Lee et al.^^ imposed numerically these symmetries in order to gain sub- 
stantial savings in computational resources. Unlike^^, our computations were performed without imposing 
any symmetry constrains, allowing thus the turbulence to evolve freely with the view that the initial TG 
vortex symmetries will break at high enough Reynolds numbers. However, even for our highest resolution 
simulations with Taylor Reynolds number of the order of 100 the TG symmetries are not broken within the 
time interval of reaching the peak of dissipation. They seem to be a strong property of the MHD equations, 
preserved by time evolution of the solutions (see alscP^. 

Due to the fact that there are special global restrictions on these TG fiows, we further consider a run 
with random initial conditions (run "R" hereafter) for comparison, ensuring that H^, = He = and kinetic 
helicity Hu = (n • u?) = at time t = 0. During the time evolution magnetic and cross helicity remain 
zero for all times relative to the total energy. However, the kinetic helicity reaches an approximate value of 
Hui/Ef < 0.2 at its absolute maximum over time but when dissipation is maximum Hui/Ef < 0.04. 

We report results based on the analysis of decaying MHD turbulence simulated with N = 1024^ grid 
points. In order to obtain the broadest inertial range, runs I, A and C are initialised at the largest scales 
and run R at wavenumbers k = 1 and 2, adding extra randomness. At time t = the fields are normalised 
such that the kinetic and magnetic energies are in equipartition, i.e. Eu{t = 0) = Ei){t = 0) = 0.125. 
Note that all fiows have unit magnetic Prandtl number (i.e. u = k). The numerical parameters of our 
computations are provided in Table [l| 

TABLE I: Numerical parameters of the DNS. The values presented are taken at the peak of total 
dissipation. Note that kmax = ^/3, using the 2/3 dealiasing rule. 



Run 


N 


V 

(xlO-^) 


Rex^ 


(xlO-^) 


At 

(xlO-^) 


(xlO-3) 


u' 


h' 


f^max'nt 


R 


1024 


5.5 


140.7 


8.33 


2.15 


7.80 


0.36 


0.48 


2.66 


I 


1024 


4.5 


121.8 


6.84 


2.03 


6.54 


0.27 


0.62 


2.23 


C 


1024 


4.5 


138.0 


6.23 


1.35 


5.82 


0.46 


0.35 


1.97 


A 


1024 


4.5 


115.1 


3.76 


1.40 


5.77 


0.37 


0.46 


1.99 



The rms velocity u' is defined as 



1/2 

u'=[^ I Eu{k)dk ) . (8) 



f- f^ 



and similarly for b', the rms of the magnetic field. The integral length scales are then defined as the total, 
kinetic and magnetic integral length scale, respectively 

_ 3n J k-^Et,u,b{k)dk 



and likewise for the Taylor scales 






In Table [ij we report the total integral and Taylor length scales as well as the Reynolds number based on 
At given by Rex^ = u'Xt/f. Finally, the smallest length scale in our flows is defined based on K41 scaling 
rjt = (y^ jti)^!^^ where e^ = i^(|u?p) + ^(|jP) is the total dissipation. The time we address in this study 
is the moment which the dissipation reaches its maximum value and therefore the highest scale separation 
occurs r]i <^ i <^ Lt^ where ^ is a typical length scale in the inertial range. Therefore, the values provided in 
Table [l| correspond to that moment. 

III. ENERGY SPECTRA 

Figure IT] presents the three-dimensional compensated total energy spectra k^Ei(]i) that we obtain at the 
peak of dissipation for all the runs of Table |lj The spectra are compensated with the scaling exponents 
p = 2, 5/3 and 3/2. The small peaks at high wavenumbers show the quality of our simulations. 
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FIG. 1: (Color online) Three-dimensional compensated total energy spectra hPEt(k) with scaling 
exponents p = 2, 5/3, 3/2 for a) run R, b) run I, c) run C and d) run A of Table llj 

According to Lee et aL^I^, the energy spectrum of the magnetically dominated fiow I, i.e. Ei)/Eu > 1 (see 
also u' and b' value s in T able S|, is close to a k~'^ power law (see Fig. flp), which is the weak turbulence (WT) 
theory expectatiorP^'^. Here, we would like to emphasize, however, that the WT scaling {Ej_ oc k^ ) is 
for an anisotropic energy spectrum, where perpendicular denotes the direction relative to an imposed large 
scale mean magnetic field ^o that does not exist in this flow. In fact, it is argued that in MHD turbulence 
there is no prescribed cascade in the parallel direction^^. This is based on the idea that small-scale turbulent 
fluctuations become anisotropic, as it is easier to shuflle strong magnetic fleld lines than to bend them due 
to the preventing action of the Lorentz force j x b. caused by a large-scale fleld Bq. 

Furthermore, Lee et al.^^ argues that the kinetic energy dominated flow C, i.e. Eij/Eu < 1 (see also 
Table yy), is compatible with a /c~^/^ slope (Fig. lib) and the less magnetically dominated flow A is near a 
^-5/3 scaling (Fig. Illi). In addition, we report that the power law of the total energy spectrum for our 
also magnetically dominated run R (see Table |l| seems to be between k~^^^ and A:~^/^. The difference 
between these two power laws is subtle enough that any type of contamination, such as intermittency or any 
dissipative small-scale effects, will blur the results. However, even a k~'^ spectrum which is slightly more 
transparent in these high Reynolds numbers can be misinterpreted. For example, in contrast to^^, we claim 



that the total energy spectrum of run A (Fig. Illi) scales like Ef oc k~^ but we leave this to the readers' 
judgement. 

Therefore, the following questions are raised: How can we circumvent this ambiguity of the results? Is 
there a dependence of small scales on the large scale initial conditions and thereby non-universality in 
decaying MHD turbulence? Since limited information can be extracted just from the slopes of the spectra, 
we try to answer these questions by examining the topology of the small scales through the invariants of 
related gradient statistics, which some have shown universal characteristics for hydrodynamic turbulent 
flows. 



IV. CLASSIFICATION OF THE FLUID FLOW TOPOLOGY 

An approach that provides a well-defined and unambiguous language to describe eddying motions and fiow 
patterns is the framework of critical point concepts from bifurcation theor}/^^, which was stud ied ex tensively 
in the context of hydrodynamic turbulent fiows by Perry, Chong, Cantwell and co-worker j-*^^*^^^ ^. Here 
we provide a brief outline on the background material related to the geometric invariants of second-order 
tensorial quantities in turbulence before going to consider various statistics of these invariants in the following 
sections. Extensive reviews on the subject can be found ini^i^^l^ and references therein. 

Geometric invariants remain unchanged under the full group of rotations (i.e. rotations plus refiections)P^, 
being independent of the frame of reference. Any traceless second-order tensor M has the following char- 
acteristic polynomial 

det[M - Xil] = ^ Af + PA,2 + QAi + i? = (11) 
where A^ are the eigenvalues of M and its invariants are 

P = -tr{M) = -(Ai + A2 + A3) = (12) 

Q=l[P^- tr{M^)] = A1A2 + A2A3 + A3A1 (13) 

R=- det(M) = -A1A2A3 (14) 
The value of the discriminant for P = is 

D=^iR^ + Q^ (15) 



and provides a general classification for the solutions of the cubic equation (11) dividing the (Q,R) space 
into the following regions 

1. D > 0: 1 real & 2 complex-conjugate eigenvalues 

2. D = 0: 3 real eigenvalues of which 2 are equal 

3. D < 0: 3 real distinct eigenvalues 

which correspond to various local fiow topologies. In this study, the first invariant is P = from definition 



( 12 ) since the vector fields that we consider are solenoidal. 



V. INVARIANTS OF THE VELOCITY GRADIENT, THE STRAIN RATE AND ROTATION RATE TENSORS 

A. Joint PDFs of the velocity gradient invariants 

The velocity gradient tensor A = Vn can be decomposed into a symmetric and skew-symmetric compo- 
nent, 

A = S ^ ft = Sij - \eijkUJk (16) 

where S = ^(Vn + Vi/-^) and Q = ^(Vu — Vn-^) are the strain rate and rotation rate tensors, respectively. 



According to Eqs. ( |T3| ) and (14), the second and third invariants of A are 

QA = l[i^^-2tr{S^)] (17) 



and 



RA = -\[tr{S'') + liOiUjjSij], 



(18) 



respectively. Here we are interested in the joint probability density function (PDF) of these invariants. A 
diagram of this joint PDF called the {Ra-,Qa) invariant map is presented in Fig. [21 labelling the various 
topological classifications. 
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FIG. 2: Diagram of the {Ra^ Qa) invariant map indicating the local fiow topologies related to each zone. 



If Qa > then enstrophy u?^ dominates over tr{S'^) and vice versa if Qa < 0. For positive values of 
Ra the topologies are unstable, whereas for negative Ra the topologies are stable. Moreover, the Da = 
line (see Fig. p|, where Da — ^Ra + Qa ^^ ^^^ discriminant, divides the invariant map into two regions. 
One where Da > with one real and two complex-conjugate eigenvalues as solutions of Eq. (11) for the 
velocity gradient tensor and the other where Da < with three real distinct eigenvalues. Note that along 
the vertical Ra = axis one of the eigenvalues is zero and therefore locally the fiow topology is invariant in 
this direction. 

Now, if Qa is much greater than zero (i.e. Da > 0) then Ra ^ —\uOiUJjSij. In this case, for Ra < 
vortex stretching dominates over vortex compression, whereas for Ra > vortex compression dominates 
(see Fig. p|. On the other hand, if Qa is much less than zero and Da < then Ra ^ —^tr{S^). In this case, 
Ra > locally is related to a sheetlike structure (or unstable node/saddle/saddle topologies according to 
the terminology of Chong et al.^ whereas Ra < with a tubelike structure (or stable node/saddle/saddle 
topologies^. This will also beco me m ore transparent when we will deal later with the third invariant of 
the strain rate tensor (see section [VB ). 

In hydrodynamic turbulent fiows, ranging from atmospheric boundary layers to free shear fiows in wind 
tunnels and even simulations of compressible turbulence, there is the prominent tendency of the joint PDF 
of {Ra^ Qa) to develop an inclined teardrop shape. This shape aligns with the second and fourth quadrants, 
with a cusp lying along the Ra > 0, Da = branch (see for example Fig. 10.1 in^^) and is considered to be 
a universal feature. Therefore, there is a preference for vortex stretching and sheet like structures. In many 
visualisations of enstrophy in hydrodynamic turbulent fiows the dominant structures seem to be tubelike 

struct ures b ut between these vortex-tubes there are sheetlike structures, where most of the dissipation is 
locatecPEH. 

Before analysing the results, we would like to note that the aspect ratio of the axes of all joint PDFs, that 
are reported in this paper, has been kept the same but the abscissa and the ordinate are different to refiect 
the change in magnitude of the plotted quantities in the four fiows that we consider. In addition, the points 
near the origin correspond to low gradient values associated with the large scale motions, whereas points 
far away characterize the high- gradient small scales. All the joint PDFs were computed at the instant of 
maximum dissipation. 

In Fig. |3]we present the joint PDFs of Ra versus Qa for all the decaying MHD runs of Table |lj The most 
important outcome from the plots in Fig. [3] is that the shape of the joint PDF of Ra with Qa is not universal 
in decaying MHD turbulence and small-scales seem to depend on the large-scale initial conditions. However, 
we should be cautious here as it is not clear whether the self-preservation of the TG vortex symmetries 
during the evolution restrict the dynamics in some way. 
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FIG. 3: (Color online) Joint PDFs of the second invariant Qa and the third invariant Ra of the velocity 

gradient tensor normalised appropriately by powers of the mean enstrophy for a) run R, b) run I, c) run C 

and d) run A of Table III The line Da = ^R\ + Qa = ^ ^^ plotted for reference. 



On the other hand, it is clear that there is a modest but still present trend of the {Ra^ Qa) map to align 
along the second and fourth quadrants for our simulation with random initial conditions (Fig. Isk). It is 
noteworthy that the shape of the joint PDF is more symmetric with respect to the Ra = axis in comparison 
to hydrodynamic flows (see for example^. Run I gives a striking joint PDF (Fig pp) with a significant 
percentage of its points lying in the first quadrant and with high absolute values of Qa in comparison to the 
rest of the runs. Points of the joint PDF in the first quadrant that are far from the origin (see Figlsb) are 
associated with very low rates of kinetic energy dissipation. This suggests that the structure is likely to be 
quite long-lived. Run C seems to resemble more the teardrop shape of hydrodynamic turbulence with the 
classic narrow cusp in Ra > 0, Da = branch (Fig. ^). Finally, Fig. 3]i shows the {Ra^Qa) map of run 
A, which has a shape with features in between the random MHD and hydrodynamic turbulence. In other 
words, there is a modest tendency of the joint PDF to align with the second quadrant like in the random 
MHD run (Fig. ^) but there is a high correlation between Ra > and Qa < values forming a long cusp 
in analogy to hydrodynamic turbulent fiows. 



B. Joint PDFs of strain rate invariants 



Setting ft to zero or essentially u) to zero in Eqs. (17) and (18), we can obtain the invariants of the strain 
rate tensor, which are 
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(19) 
(20) 



The (i?s, Qs) invariant map features the geometry of the local straining of the fiuid elements (see Fig. [4|. 
The second invariant Qs is related to viscous dissipation e = 2h'S'^ through Qs = —\^/i^ because the strain 
rate tensor is symmetric, i.e. S'^ = SijSji = tr{S'^). So, locations with Qs much less than zero are highly 
dissipative regions. Note that Qs is negative definite. The third invariant Rs has two important physical 



Tube structure 
Rs<0 

Axisym. 
contract 
Ds=0 

\ / 
|2:-1:-1 / Hydro 

turbul. ^^ 

2D flow 



Sheet structure 
Rs>0 

Axisym. 
expansion 



1:0: 





R. 



FIG. 4: Diagram of the (i?s, Qs) invariant map. Each plotted curve corresponds to the fohowing flow 

geometries: Ai : A2 : A3 = 2 : — 1 : — 1 (axisymmetric contraction), 1 : : — 1 (two-dimensional flow), 

3:1: —4 (biaxial stretching) and 1 : 1 : — 2 (axisymmetric stretching). 



meanings. First, it is proportional to strain skewness SijSjkSki, which appears as part of the production 
term in the evolution equation of S'^ (see^^). Second, it can be written as a function of the eigenvalues of 



Rg 



3(^1 + ^2 + ^3) — ~^1^2A3 



(21) 



since tr{S) = Ai + A2 + A3 = due to incompressibility, with Ai > A2 > A3. Owing to the symmetry of 
Sij all eigenvalues are real and therefore the {Rs^Qs) invariant map is contained only in the region where 
Dg = ^R^-\-Ql < (see Figs. Hand^. So, Rs > implies production of S'^ and hence of viscous dissipation 
with Ai, A2 > and A3 < related to sheetlike structures. On the contrary, Rs < indicates destruction of 
S^ with Ai > and A2, A3 < associated with tubelike structures. Note, therefore that sgn{Rs] 
We should point out here that if we deflne the following ratio a = A2/A1 of the eigenvalues of Sij, then each 



sgn(A2 



value of a corresponds to a line in the (Rs^Qs) plane with the following expression 



Rs 



{-Qsf^^a{l + a){l + a + a^)-^/^ 



(22) 



where each line is associated with a flow topology (see caption of Fig. EP^J^. 

The {Rs^Qs) invariant map in many hydrodynamic turbulent flows away from boundaries manifests a 
tendency for the Rs > region, implying a predominance of sheetlike structures related to the strain 
rate (see for example Fig. 8c iiP^. In particular, numerical and experimental evidences in homogeneous 



3:1:- 



hydrodynamic turbulence propose the ratios of the mean eigenvalues of Sij to be (Ai):(A2):(A3) 
(see the corresponding line in Fig. [4|. 

The joint PDFs of Rs versus Qs for the four runs of Table |l] are illustrated in Fig. [5] Their dependence on 
initial conditions is clearly depicted. The shape of the {Rs^Qs) iiiap for run R (Fig. Isk) moves away from 
the Ds = line towards the Rs = axis expressing a more quasi two-dimensional (2D) character of the 



structures related to Sij than in hydrodynamic turbulent flows away from the boundaries. We should point 
out, however, that this particular shape is reminiscent to the joint PDFs of {Rs^Qs) found in the buffer 
layer, i.e. a region very close to the wall, of wall-bounded turbulent shear flows (see for example Fig. 6f 
iii^. The joint PDF of run I (Fig. pp) is aligned along the Rs = with some highly dissipative small scales 
in contrast to the rest of the runs. The local topology in this case seems to have a strong tendency towards 
quasi two-dimensionality. Part of the shape of this joint PDF can be explained through two-dimensional 
shearing (or vortex sheet), i.e. 



^0 dyU:, 0^ 




.0 a. 



(23) 



0; 



which gives Qs = —\[{dyUx)'^ + i^yUz)'^] and i^^ = in analogy to the influence of the wall on the velocity 
gradient in wall-bounded flows^. The {Rs, Qs) invariant map of run C in Fig. pp also falls away from the 
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FIG. 5: (Color online) Joint PDFs of the second invariant Qg and third invariant Rs of the strain rate 

tensor normalised appropriately by powers of the mean enstrophy for a) run R, b) run I, c) run C and d) 

run A of Table |lj The line Ds = ^R^ ^ Q^ = is plotted for reference. 



Dg = branch with low correlations between Rs and Qg. Finally, the joint PDF of run A (Fig. 
almost identical in shape but less correlated with respect to Fig. [5k. 
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We now try to summarise and clarify our arguments by tabulating the mean eigenvalues of the strain rate 
tensor and their ratios for all our runs in Table [ll] but also by plotting the curves that can be constructed 
from Eq. (22) using the mean eigenvalues of Table [TT| (see Fig. IgI). 



TABLE II: Mean eigenvalues of the strain rate tensor Sij and their ratios. 



Run 


(Ai) 


(A2) 


(A3) 


(Ai) : (A2) : (A3) 


R 


0.14 


0.03 


-0.17 


5 


1 


-6 


I 


0.25 


-0.00 


-0.25 


1 





-1 


C 


0.25 


0.04 


-0.29 


6 


1 


-7 


A 


0.25 


0.04 


-0.29 


6 


1 


-7 



In Fig. [6) we plot for reference the curve that corresponds to 3:l:-4, the characteristic eigenvalue ratios for 
homogeneous hydrodynamic turbulent flows that we denote as "HD". In that respect, all the ratios of the 
mean eigenvalues that we obtain are different than 3:l:-4. However, all the cases represent biaxial expansion 
apart from run I, which is characterised by quasi two-dimensionality with weak biaxial contraction (see 
Table In] and Fig. l6|. Figure Ig] makes clear that on average the flow topologies related to Sij of run C and A 
are close to run R giving weight to our argument for the similarity of their (i?s, Qg) joint PDFs. The curve 
for run I also summarises Fig. [SJd by demonstrating that the quasi 2D structures associated to the strain 
rate tensor are weakly contracted in a average sense. 
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FIG. 6: (Color online) Plots of Eq. (22) using the mean eigenvalues of Sij from Table pi) The dashed line 



Ds = '4r' 



^^ = is plotted for reference. With HD we label the curve that corresponds to 3:l:-4, the 
characteristic eigenvalue ratios for homogeneous hydrodynamic turbulence. 



C. Joint PDFs of the second invariants of the strain and rotation rate tensors 

Another important joint PDF to analyse is the one of —Qs versus the second invariant of the rotation 
rate tensor, Qa;, which is in fact the only invariant for fl. To see this, set S to zero in Eqs. (17) and (18), 
then 



= -hr{n') = i 



CJ 



(24) 



which is positive definite and it is related to the second invariants of A and S through Q^ = Qa — Qs- 
The (Qu;, —Qs) invariant map that is shown schematically in Fig. ItI identifies the relative importance of the 
straining and rotational part of velocity gradient tensor. A good example that describes simply the physical 




FIG. 7: Diagram of the (Qc^, —Qs) invariant map pointing out the important regions related to strain and 

rotation. 

meanings of Fig. [7|is the Burger's vortex tube^. As it was mentioned before Qs characterises the topology 
associated with viscous dissipation. So, points near the — Qs axis refiect nearly pure straining motions, i.e. 
regions of strong dissipation but negligible enstrophy, like outside and away from the Burger's vortex tube. 
On the other hand, points close to the Q^j axis are in nearly pure solid-body rotation, like at the centre of 
the Burger's vortex tube with high enstrophy but very weak dissipation. Regions with comparable strain 
rate and rotation map to points close to the Q^ = —Qs line, which correspond to vortex sheets. 

Generally, from observations in many hydrodynamic turbulent fiows, regions of intense enstrophy tend to 
be concentrated in tubelike structures, whereas regions of high dissipation are not correlated with regions 
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of concentrated enstroph};^^. So, the joint P DF of Q^j versus —Qs is very spread for many hydrodynamic 
turbulent flows away from walls (see results irr^^). 

Figure [8| shows the joint PDFs of Q^ versus — Qs, normalised with the mean enstrophy, for the four runs 
of Table lir The dependence on initial conditions is pronounced once more in these plots. The {Qoj-,—Qs) 




Q /<|co|2> 



(d) Run A 



FIG. 8: (Color online) Joint PDFs of the second invariants of strain rate and rotation rate tensors 
normalised by the mean enstrophy for a) run R, b) run I, c) run C and d) run A of Table llj 

invariant map of run R (Fig. ^) is very different to hydrodynamic turbulence away from walls. Here, the 
joint PDF is concentrated around the Q^ = —Qs line demonstrating stronger correlation between these 
two variables. This result in conjunction with the outcome from Fig. [6] conflrms many visualisations of 
homogeneous MHD turbulent flows^^, which illustrate large population of sheetlike rather than tubelike 
structures. 

The shape of the joint PDF (Qu;, ~Qs) for run I is even more extreme with a very narrow distribution 
along the main diagonal (Fig. jsj)), where regions of high dissipation are strongly correlated by high levels of 
enstrophy particularly for points far from the origin. The high gradients in this flow can be well approximated 
by Eq. (23) where Q^ = —Qs = \[{dyUx)'^ + {dyUz)'^]. According to Cantwell^^ the presence or absence 
of points very far from the origin, associated with quite long-lived structures, is closely related to the 
regularity of the initial conditions. He further mentions that such structures are much less prominent in 
a flow with randomised initial conditions. Here, this is transparent if one compares the run with random 
initial conditions (Fig. |8^) with run I (Fig. JsJd). Moreover, it could be argued that the core of the joint 
PDF (Quj , ~Qs ) of run I is similar to the joint PDF obtained in the buffer layer of wall-bounded flows (see 
results b}E^^^. The {Qojj—Qs) i^ap of run A (Fig. ^) resembles Fig. ^ but with weaker correlations 
between high dissipation and high enstrophy regions. Finally, the joint PDF of Fig. (Sj), which corresponds 
to run C, presents the weakest correlations between Q^ and —Qs among the four cases with a weak trend 
of alignment along the main diagonal. 



D. Flow structures and enstrophy dynamics 



Various flow fleld quantities were viewed interactively using a visualisations software^ to get an idea of 
the spatial structures in our flows. In order to substantiate our approach, we present indicatively plots of 
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iso-contours of the vorticity field in our [0,27r]^ periodic boxes at the moment of maximum dissipation for 
the four runs of Table [l| (see Fig. |9|. Figure |9|i (run R) displays iso-contours of vorticity for |a;| > 3a;' 




(a) Run R 



(b) Run I 




(c) Run C 



(d) Run A 



FIG. 9: (Color online) Vorticity field iso-contours with \uj\ > 3a;' for a) run R, b) run I, c) run C and d) 

run A of Table |lj 

where a;' = (|u;p)^/^. The predominant structures in this plot are randomly oriented sheetlike structures 
in support of our joint PDF analysis. In comparison to the randomly oriented structures of run R, the TG 
vortex symmetries become apparent in Figs. [9]3, c and d revealing their preservation in time. Remember 
that we did not impose any symmetries during the evolution of our runs. According to the above analysis, 
the peculiar run I should be prevailed by quasi two-dimensional sheetlike structures, which are shown in Fig. 
[9p. These flat structures are formed on the insulating faces of the [0,7r]^ boxes and on their mid-planes in 
the vertical direction, i.e. z = 7r/2. The structures of run A (Fig. |9]i) are also sheetlike but more randomly 
oriented in contrast to run I. In the end, run C is a more complicated TG flow. This is demonstrated in 
Fig. [9J3 for I a; I > 3a;'. It is interesting that the initial conditions of the TG velocity with the TG magnetic 
fields for the insulating runs I and A create less randomness in the flow fields, which are mainly dominated 
by quasi 2D sheetlike structures in contrast to run C. 

According to Jimenez et al.^^, in hydrodynamic turbulent flows away from walls, it is qualitative clear 
that there is no other way of production of enstrophy other than straining of weak vorticity to form stronger 
vortex regions. Then, strain itself is induced by vorticity and the process may become non-linear. This 
mechanism is called self-amplification of velocity derivativea^^'^. 

In order to have an initial picture of this mechanism and in particular of the formation of the vorticity 
fields in our MHD flows, we examine the rate of vortex stretching 



S = 



uj ■ S ■ uj Rs — Ra 



u; 



Qu 



(25) 



which is essentially the part of the strain that is aligned with the local vorticity and it is the term that 
stretches or compresses the vortex lines in the evolution equation of the enstrophy 



dti^iv'^) =uj - S - uj ^vuj - Acj + cj • V X (j X 6). 



(26) 
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Notice that S can be written as a function of the invariants Ra^ Rs and Q^j (see Eq. (25)). 

Figure [To] shows joint PDFs of essentiahy the enstrophy (i.e. Q^) with the rate of vortex stretching S 
appropriately normahsed for ah the flows of Table [Tj Various common features can be observed in Fig. 10 




Z/<|0)|2>1/2 

(a) Run R 



Z/<|C0|2>1/2 

(b) Run I 




Z/<|0)|2>1/2 

(c) Run C 



Z/<|0)|2>1/2 

(d) Run A 



FIG. 10: (Color online) Joint PDFs of the second invariant of the rotation rate tensor Q^ and the vortex 
stretching rate S normalised appropriately by powers of the mean enstrophy for a) run R, b) run I, c) run 

C and d) run A of Table llj 

To be more specific, the highest values of enstrophy are associated with positive but low values of S, i.e. 
stretching of vorticity, whereas high rates of stretching as well as compression correlate with regions of low 
Q^. So, there is little evidence of self-stretching by structures in the fiow which have large enstrophy in 
analogy to hydrodynamic turbulence^^ ^^. Another common feature in all the plots of Fig. 10 is the tilt 
towards positive values, i.e. vorticity vectors are being more stretched than compressed. 

On the other hand, quantitative differences are evident, such as the asymmetry of the (S, Q^) joint PDFs, 
which seems to be different for each fiow. In other words, the joint PDF of run C (Fig. p!Op) is shifted more 



towards S > values, akin to hydrodynamic turbulence (see results irP^Eo^), in comparison to run A (Fig. 
[lOJi) which is closer to the joint PDF of run R (Fig. 10 i). Another quantitative difference between the four 
fiows is the very high values of enstrophy {Q^ :^ 6(|u5f)) that are obtained in run I (Fig. 10 3) for values of 
vortex stretching rate of the same order for all the fiows (i.e. S < 0.3(|u;p)^/^). 

Another important mechanism for amplification or reduction of enstrophy that exists only in MHD tur- 
bulence is that due to the Lorentz force. This process essentially manifests from the last term of Eq. (26), 
which we write here as 



uj -V X {j X b) 



(27) 



SO that it is comparable with S (see Eq. ([25|). In order to shed light on the dynamics of this term with 
respect to the enstrophy, we consider in Fig. |ll|the joint PDFs between L and Q^, normalised appropriately, 
for the four runs of Table [H 



It is characteristic for all the plots of Fig. 11 that there is a preference for L > for most of the local 
topology in the fiow. It is also common in all the four cases that the highest values of Q^ are associated 
with regions of low but positive L, whereas high values of |L| are related to regions of low enstrophy in a 
similar fashion to the self- amplification mechanism. Once more, the quantitative differences between the 
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L/<|a)|2>i/2 
(a) Run R 



2 

L/<|0)|2>1/2 

(b) Run I 




L/<|0)|2>1/2 

(c) Run C 



L/<|0)|2>1/2 

(d) Run A 



FIG. 11: (Color online) Joint PDFs of the second invariant of the rotation rate tensor Q^ and L normalised 
appropriately by powers of the mean enstrophy for a) run R, b) run I, c) run C and d) run A of Table llj 



plots of Fig. 11 are evident with the most notable being the joint PDF of run I (Fig. flTb) with the highest 
values of Q^ in terms of L. 

A comparison between the two mechanisms of amplification and reduction of enstrophy reveals the cause 
of high and low enstrophy regions. For the MHD flow with random initial conditions, L is more correlated 
with regions of higher enstrophy than S but the opposite is true for the TG flows. On the other hand, the 
lowest enstrophy regions present correlations with higher absolute values of L than S for all the runs. 



VI. INVARIANTS OF THE MAGNETIC FIELD GRADIENT, THE MAGNETIC STRAIN RATE AND THE 
CURRENT DENSITY RATE TENSORS 

In this section we try to classify the topology related to the magnetic field by extending the above joint 
PDF analysis for the invariants of magnetic field gradient tensor as well as for the invariants of its symmetric 
and skew- symmetric components. 



A. Joint PDFs of the magnetic field gradient invariants 



The magnetic field gradient X 
component 



V6 can be also decomposed into its symmetric and skew-symmetric 



X = K + J = K^p 



)^a(3^J'y 



(28) 



where K = ^(V6 + V6-^) and J = ^(V6 — V6-^) are the magnetic strain rate and current density rate 
tensors, respectively. The skew-symmetric part of X is related to the electric current through Ampere's law 
V X 6 = fioj where /io = (^cr)~^ is the permeability of free-space and a is the electrical conductivity. When 
the magnetic field lines are bended, current is produced providing a Lorentz force that inhibits the bending 
of the field lines. On the other hand, the symmetric part of X characterises the force-free regions in the 
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magnetic field, where j = and therefore j x b = 0. An important relation one can easily derive by taking 
the divergence of Eq. ^ and using the fact that our fields u and b are solenoidal is the following Poisson 
equation 



V^P 



V-[( 



U X (jj) 



(i X b)] 



— \^a(3 ~ ^a(3) + \^a(3 ~ '^afdJ' 



(29) 



What is interesting in this expression is the interchange between the symmetric and skew-symmetric tensors 
of Vn and V6 related to V^P. It is also appealing that the viscous dissipation is related to the symmetric 
part of the velocity gradient, whereas the Ohmic dissipation to the skew- symmetric part of the magnetic 
field gradient. 

Now, we consider the joint PDF of the second and third invariants of X, which are defined according to 
Eqs. (pi) and (pi) as follows 



Q 



X 



Rx - -3 
For the classification of the magnetic field structures the Dx 



l[f-2tr{K')] 


(30) 


r{K') + IjaJpK^p]- 


(31) 


ictures the Dx = ^Rx + Qx 


= line was included in 
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FIG. 12: (Color online) Joint PDFs of the second invariant Qx and the third invariant Rx of the magnetic 
field gradient tensor normalised appropriately by powers of the mean squared current density for a) run R, 



b) run I, c) run C and d) run A of Table I The line Dx 






27 E)2 
-4^X 



Qx = is plotted for reference. 



the plots of Fig. [12] The topological classification emerging from the joint PDFs of Rx and Qx can be 
interpreted in analogy to the invariant map of the velocity gradient (Fig. [2|. Note, however, that the 
individual terms of the third invariant in Eq. (31) do not appear in any evolution equation. Thus, Rx does 



not have a physical meaning here but it is mathematically important for the classification of the magnetic 
field structures, in terms of the eigenvalues of Xaf3 associated with these structures. 

In contrast to the invariants of the velocity gradient, the {Rx^Qx) invariant map does not show a 
particular tendency towards any quadrant (see Fig. 12). For all the runs the core shape of the joint PDF is 



16 



symmetric along the Rx = axes, meaning that there is a balance between stable and unstable structures. 
The small scales, on the other hand, are slightly different especially for run I (Fig. ^^) and run C (Fig. 
T2fc). Moreover, the joint PDF for runs C and A (see Fig. 12 3 and d, respectively) diminish towards the 



origin of the axes. In general, one could claim that this symmetric shape seems to be a general characteristic 
for the magnetic field gradient for all initial conditions with some small deviations, which might be due to 
the TG vortex symmetries. 



B. Joint PDFs of the magnetic strain rate invariants 



Looking at the joint PDFs of the second and third invariants of K we can study the geometry of the local 
magnetic straining. The invariants of the magnetic strain rate tensor can be obtained by setting j = in 
Eqs. (30) and (31), which reduce to 



QK = -^tr{K^) 



(32) 



and 



Rk = -\tr{K^), 



(33) 



where Qk is negative definite due to the symmetric nature of K. Note that Qk is not directly related to 
Ohmic dissipation in contrast to the Qs for viscous dissipation. Then, the physical interpretation of the 
(RkjQk) invariant map is quite different from Fig. [4] but similar in terms of flow topology. So, very low 
values of Qk in Fig. 13 can be physically interpreted as regions of high magnetic-strain or regions where the 
Lorentz force is small. The third invariant Rk can be written as the product of the eigenvalues of K^is in 
analogy to Rs (see Eq. (21)). Then, the interpretation of Rk in terms of sheetlike and tubelike structures 
is also determined by sgri(!Rx) = sgn(A2). 

The joint PDFs between Rk and Qk^ representing the local topology of the structures related to magnetic 
strain rate, appear to be symmetric along the Rk = axis for most of the runs of Table |l] (see Fig. 13). 



In detail, the joint PDF of run R (Fig. 13 1) illustrates an equipartition between tubelike and sheetlike 



structures associated with K. The shapes of the [RkiQk) invariant map for runs I and A (see Figs. 13 3 
and[l3]i respectively) are also symmetric and they can be well approximated by a magnetic field gradient of 
the form of Eq. (23) with Qk = —\[{9ybx)'^ + i^ybz)'^] and Rk = 0. In Fig. 13 3 there are very low values 
of Qk correlated with Rk = in comparison to the rest of the flows. Therefore, this approximation for 
X is especially valid for the small scale structures that correspond to low values of Qk in this joint PDF. 
Figure fl3fc (run C), on the other hand, is slightly asymmetric, showing a tangible inclination of the joint 
PDF towards Rk < 0. This implies that there is a preference for the intermediate eigenvalue of Kai3 to be 
negative and hence a tendency for tubelike structures. 

Now, we attempt to provide an outline of the joint PDFs of Fig. [l3]by tabulating the mean eigenvalues 
of the magnetic strain rate tensor (see Table III ) and by plotting the analogous expression to Eq. ( 22 ) for 
Rk and Qk using the mean eigenvalues of Table III (see Fig. 14). The values of the mean eigenvalue 



TABLE III: Mean eigenvalues of the magnetic strain rate tensor Kai3 and their ratios. 



Run 



(Ai> 



(A2> 



<A3> 



(Ai) = (A2) = (Xs) 



R 


0.26 


0.00 


-0.26 


I 


0.44 


0.00 


-0.44 


C 


0.35 


-0.02 


-0.33 


A 


0.44 


0.00 


-0.44 



1 


: 


-1 


1 


: 


-1 


18 : 


-1 


-1 


1 


: 


-1 
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ratios tell us that on average run I and C are described by quasi two-dimensional structures in agreement 
with the joint PDF analysis. Moreover, the values 1:0:-1 that we obtain for run R agree with the argument 
that the joint PDF of Fig. [l3^ is symmetric but also express that in an average sense the flow topology is 
locally invariant in one direction. The only case that deviates from two-dimensionality is run C, which is 
on average characterised by biaxial contraction (i.e. (A2) < 0) and thereby tubelike structures. 
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(a) Run R 



(b) Run I 




(c) Run C 



(d) Run A 



FIG. 13: (Color online) Joint PDFs of the second invariant Qk and the third invariant Rk of the 
magnetic strain rate tensor normalised appropriately by powers of the mean squared current density for a) 
run R, b) run I, c) run C and d) run A of Table pj The line Dk = x^i: ^ ^k = is plotted for reference. 
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FIG. 14: (Color online) Plots of Eq. (22) using the mean eigenvalues of Ka(3 from Table III The dashed 



line Dk = x^ir ^ ^x = is plotted for reference. 



C. Joint PDFs of the second invariants of the magnetic strain and current density tensors 



The skew-symmetric part of the magnetic field gradient tensor, J has only one invariant in analogy to 
the rotation rate tensor ft. This can be obtained by letting K to be zero in Eqs. (30) and (31), then 



which is also related to the second invariants of X and K through Qj = Qx — Qk- 



(34) 
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The (Qj, —Qk) invariant map describes the relative importance between the straining and rotational part 
of the magnetic field gradient in analogy to (Qa;, —Qs) map for the velocity gradient (see Fig. [t]). However, 
the important difference in this case is that the rotational part of X is directly related to Ohmic dissipation 
and not the straining part. Hence, the points of the joint PDFs close to the Qj axis that are nearly in 
solid-body rotation are regions in the fiow of strong Ohmic dissipation and negligible magnetic straining in 
contrast to the picture we get from Fig. [7| On the other side, points adjacent to the —Qk axis express 
nearly pure magnetic straining motions in regions of where the current is negligible and thereby Lorentz 
force is suppressed. 



The joint PDFs of Fig. 15 show that points near the axes are rare in MHD turbulent fiows and are related 



only to the large scales of the flows, where Qj and —Qk are small in comparison to (|j| ). Most of the 
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FIG. 15: (Color online) Joint PDFs of the second invariants of magnetic strain rate and current density 
rate tensors normalised by the mean squared current density for a) run R, b) run I, c) run C and d) run A 

of Table [Tj 



points in the plots of Fig. 15 lie near the main diagonal, revealing that Ohmic dissipation occurs in current 
sheets. Here, the magnetic field gradient tensor can be well approximated by the form of Eq. (23), which 
gives Qj = —Qk = \[{9yhx)'^ + i.^y^z)'^]- This is particularly a good approximation for runs I and A (see 



Figs. 15 3 and 151 respectively), where Qj and —Qk are strongly correlated for all scales. It can also be 



argued that this approximation is also valid for the small scales of runs R and C (see Figs. 15ti and 
respectively) that correspond to high values of Qj and -Qk- 
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D. Structures In the current density field 



To further validate our joint PDF approach, we present indicatively iso-contours of current density (Fig. 
16) in our [0,27r]^ periodic boxes at the moment of maximum dissipation for all the runs that we have 
considered (see Table ^. All the visualisations of Fig. 16 display current density iso-contours with \j\ > 6j' 
where / = (|jp)^/^. The field of current density for run R (Fig. 16 i) is composed by randomly oriented 
sheetlike structures which seem to be extremely thin, supporting the fact that the values of the mean 
eigenvalue ratios for the magnetic strain rate tensor are 1 : : — 1. It is clear that locally quasi two- 
dimensional structures are the dominant structures in Figs. 16 3 and [l6]i (runs I and A respectively). 



19 




(a) Run R 



(b) Run I 




(c) Run C 



(d) Run A 



FIG. 16: (Color online) Current density field iso-contours withjjl > 6/ for a) run R, b) run I, c) run C 

and d) run A of Table III 



validating the joint PDFs of {Rr^ Qk)- These 2D current sheets are also the structures where most of the 
Ohmic dissipation occurs. For run I these dominant structures are formed at the faces of the [0,7r]^ boxes, 
whereas for run A these are randomly oriented. On the other hand, run C (Fig. 16 3) seem to be dominated 



by tubelike structures but one can also observe the coexistence of isolated thin current sheets in agreement 
to our analysis. Finally, the TG vortex symmetries are clearly depicted in these visualisations with each TG 
flow having different degree of randomness. This raises again questions as to what degree these symmetries 
restrict the dynamics of the flows. 

Here, we would like to emphasize that the structures related to the magnetic fleld gradient have different 
characteristics than those related to the velocity gradient. This might well be a reason that the energy 
spectra that we obtain, as well as other studies, for the kinetic and magnetic energy (not shown here) seem 
to obey different scaling exponents. Moreover, we conjecture that these quasi two-dimensional organised 
structures that appear in runs I and A both in the vorticity and the current density flelds are the reason to 
obtain a k~'^ scaling that we observe in the total energy spectra in Figs, [ip and [ip. This, however, needs 
to be further investigated and it will be reported elsewhere. 



VII. CONCLUSIONS 



The universality of the energy spectrum in MHD turbulence is in doubt by various studies. One aspect 
is the manifestation of different, dubious scaling exponents. In order to avoid ambiguity between scaling 
exponents, we explore various statistics based on the invariants of the velocity gradient and related tensors. 
Note that for a big family of hydrodynamic turbulent flows, the joint PDF of the invariants of the velocity 
gradient is generally considered to be universal. We further extend this analysis to the invariants of gradient 
statistics related to the magnetic fleld. In particular, we explore DNS data of decaying MHD turbulence 
with random initial conditions as well as a set of three different Taylor- Green type initial conditions without 
imposing any symmetry constrains in our flows during their evolution. The TG flows were chosen to 
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be examined since recently, Lee et al.^^ reported that the scahng of the energy spectrum at the peak of 
dissipation depends on the initial conditions. 

Our study attempts to classify the structures of our MHD flows. The structures related to the strain rate 
tensor are predominantly sheetlike structures (i.e. (A2) > 0) for all the flows apart from run I (see Fig. [6|, 
which is quasi two-dimensional (i.e. (A2) — 0). The biaxial stretching for our MHD flows is different in 
comparison to hydrodynamic turbulence, namely (Ai) : (A2) : (A3) = 3 : 1 : — 4 (see Table [Tl[) . Furthermore, 
the enstrophy dominated regions are well correlated with region of high viscous dissipation in contrast to 
hydrodynamic flows. We also find that viscous dissipation is an intrinsic element of vortex sheets. 

On the other side, magnetic field consists of quasi two-dimensional structures, i.e. (A2) = 0, for all the 
cases apart from run C, which is on average dominated by tubelike structures, i.e. (A2) < (see Table 



III). The correlation between magnetic strain dominated regions and regions of high Ohmic dissipation is 
generally stronger than the correlation between enstrophy and viscous dissipation. We also obtain that 
Ohmic dissipation resides in current sheets, which are thinner than the vortex sheets in the same flow. 
Visualisations support further our joint PDFs analysis of the invariants. 

Our results also demonstrate that small scales depend on the initial conditions in decaying MHD tur- 
bulence. This is dramatically illustrated through the joint PDF of Ra with Qa (see Fig. |3|, which has 
a universal teardrop shape for hydrodynamic turbulence away from walls. Lack of small scale universality 
in decaying MHD turbulence will have important implications in modelling. The main idea of small-scale 
universality applied in LES (i.e. that although large scales may be dependent on boundary conditions or 
initial conditions, smaller scales are less flow dependent and more amenable to modelling) seems to fail for 
decaying MHD turbulent flows. Therefore, if MHD turbulence is non-universal, then the construction of 
subgrid-scale models for MHD flows might be doubtful. 

However, there is an important element regarding the TG flows that one has to address before claiming 
that small scale universality is absent in these flows. This element is the self-preservation of TG vortex 
symmetries during the evolution of the flow, which seem to be a strong property of the MHD equations. 
Therefore, a natural question that emerges is: what happens if we perturb the TG flows in order to break 
these symmetries before the peak of dissipation? Will the joint PDFs converge to a single shape or/and 
the scaling of the energy spectra to a single value? What is the role of the symmetries imposed by the 
initial conditions in terms of the dynamics? Do we have classes of universality for these moderate Reynolds 
numbers or is there a universal power law in the high Reynolds number limit? We plan to address these 
questions in our future work. 
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